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Abstract 


Modern multi-layer insulation (MLI) allows to sharply reduce the heat leak into cryogenic propellant 
storage tanks through the tank surface and, as a consequence, significantly extend the storage 
duration. In this situation the MLI penetrations, such as support struts, feed lines, etc., become one 
of the most significant challenges of the tank’s heat management. This problem is especially acute 
for liquid hydrogen (LH2) storage, since currently no efficient cryocoolers exist that operate at very 
low LH2 temperatures (^20K). Even small heat leaks under microgravity conditions and over the 
period of many months give rise to a complex slowly-developing, large-scale spatiotemporal physical 
phenomena in a multi-phase liquid- vapor mixture. These phenomena are not well-understood nor 
can be easily controlled. They can be of a potentially hazardous nature for long-term on-orbital 
cryogenic storage, propellant loading, tank chilldown, engine restart, and other in-space cryogenic 
fluid management operations. To support the engineering design solutions that would mitigate these 
effects a detailed physics-based analysis of heat transfer, vapor bubble formation, growth, motion, 
coalescence and collapse is required in the presence of stirring jets of different configurations and 
passive cooling devices such as MLI, thermodynamic vent system, and vapor-cooled shield. To 
develop physics-based models and correlations reliable for microgravity conditions and long-time 
scales there is a need for new fundamental data to be collected from on-orbit cryogenic storage 
experiments. Our report discusses some of these physical phenomena and the design requirements 
and future studies necessary for their mitigation. Special attention is payed to the phenomena 
occurring near MLI penetrations. 
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1 Introduction 

This paper aims at discussing some of the basic physics issues associated with long-term storage 
of cryogenic liquids in zero gravity or microgravity environments. By “long-term” we mean, for 
example, the durations of the currently envisioned extended storage periods in the low earth orbit 
(LEO), which range from months to years. For NASA’s present and future space exploration 
missions, understanding the behavior of cryogenic liquids over long periods of storage is of crucial 
importance, because of the fundamental role played by cryogenic propellants, primarily liquid 
hydrogen (LH2) and liquid oxygen (LOX), in rocket propulsion, specifically for long-range missions 
[1]. The very feasibility of using liquid propellant engines based on LH2 and LOX in the future 
long-range missions depends on the success of storing these propellants under microgravity or zero 
gravity for extended periods of time. One of the currently considered exploration strategies calls 
for the development of propellant storage and transfer facilities in LEO [2]. These “fuel depots” 
will need to be able to spend significant amounts of time (at least on the order of several months) in 
LEO without any substantial propellant losses due to boil-off [2-5] . With the current passive heat 
insulation technologies, it is theoretically possible to reduce the cryogen boil-off rate to below 3% 
per month [6]. Even so, this issue becomes a challenge when the required storage duration exceeds 
6 months, and yet a greater one for manned missions to Mars [7]. 

Cryogenic fluid management (CFM) in microgravity provides a number of fundamental physical 
challenges, many of which were previously discussed in the literature [6,8-13]. This is especially 
relevant to storage of LH2 because of its low critical temperature. One of the main features of 
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microgravity environments is that due to much-reduced levels of ^-forces and their generally time- 
varying character, the vapor bubbles that form as a result of boil-off at the tank walls mainly near 
MLI penetrations (hot-spots) will not rise quickly towards the ullage space, as they do under normal 
gravity. Instead, they may slowly grow to very large sizes (tens of centimeters), or they may detach 
from the wall, migrate toward the stagnation areas of the stirring flow and accumulate there, forming 
regions of saturated liquid and complicated foam-like vapor-liquid structures whose properties may 
be not easy to control. These processes are governed by the complex heat transfer mechanisms in 
the near-wall region; capillary and ^-forces; complex dynamics of nucleate boiling; bubble growth, 
detachment and collapse; chemical traces that can accumulate in the liquid with time and affect 
its properties. In long-term storage missions foam or bubble colonies can grow at the expense of 
the single ullage space. They may not be easily removed by tank pressurization because the heat 
released from vapor condensation may raise the temperature of the liquid surrounding the bubbles 
to the saturation temperature at the higher pressure. Capillary forces may be sufficiently strong 
to prevent the detachment of the foam from the tank walls by any realistic stirring flow that keeps 
the ullage intact. Basic challenges, therefore, include control of the tank pressure, temperature, 
ullage space size and location, boil-off venting, and work of liquid acquisition devices (LAD) that 
can be clogged by the foam. Similarly, since the role of buoyancy-driven convection, which is 
the main mechanism of heat transfer on earth, is greatly reduced in microgravity, heat transfer 
mechanisms will be significantly altered. Vapor and fluid motion, in turn, will be dominated by the 
capillary forces, heat transfer-mediated bubble dynamics, bubble coalescence, Ostwald ripening and 
the induced thermocapillary convection. The resulting bubble patterns and near-wall dynamics, 
especially around the MLI penetrations can substantially depend on the type of the wall material, 
chemical traces, vibrations and other external factors. 

In view of these complications, the basic technical issues that need to be dealt with in today’s 
design of successful cryogenic storage and transfer devices for long-term operation in microgravity: 
heat transfer management, pressure control, design of tank stirring, mass gauging, liquid acquisi- 
tion, and fluid transfer are much more challenging [11], compared to the Apollo era short duration 
missions, in which a low level of gravity was propulsively maintained [14-16]. 

A successful treatment of the pressing technical issues of cryogen management in microgravity 
is impossible without a thorough mechanistic understanding of the underlying physical processes of 
nucleate boiling. Surprisingly, detailed physical understanding of nucleate boiling phenomena is still 
lacking today (see e.g. [17]). This may be due to the fact that boiling is a strongly non-equilibrium 
phenomenon in which an interplay between stochastic nucleation events at the micro-scale and 
complicated deterministic nonlinear dynamics at macro-scale takes place (for reviews, see [18-21]). 
At the same time, in view of the fact that microgravity presents quite a different environment 
compared to the usual environment on earth, one should exercise caution in applying the engi- 
neering correlations developed under earth gravity conditions to the design of cryogenic systems 
to be operated in space [21,22]. To address the above technology gaps, it is necessary to collect 
fundamental data on liquid- vapor structure and dynamics during long-term storage in microgravity 
from carefully designed in-space long-duration experiments. Experimental work should be done in 
combination with a detailed physics analysis, mechanistic modeling, first principles computational 
and multi-scale approaches. 

Here we perform some basic physical estimates in order to evaluate the relative importance of 
different physical processes during long-term cryogenic storage. We concentrate our efforts on LH2, 
since it is the cryogen of primary importance to rocket propulsion and is also the most difficult 
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Figure 1. The schematics of a cryogenic storage tank (from [6]). 


in terms of CFM due to its low boiling point. Let us emphasize that we aim at obtaining only 
relatively rough estimates that take into account the long-term nature of storage. Thus, our main 
tool will be dimensional analysis, with minimal reference to more advanced mathematical tools. 
Once the main physical processes acting on the considered long timescales are identified, relevant 
space experiments can be designed and more precise calculations may be made using advanced 
mathematical and high-fidelity computational tools. In short, our main goal is to identify these 
processes and the issues, such as safety hazards and design optimization parameters, which arise 
specifically during extended periods in microgravity. 

2 Background 

We start with some basic considerations relevant to large-scale cryogenic storage in microgravity. 
A conceptual representation of a cryogenic fuel tank [6] is shown in Fig. 1. To fix ideas, let us 
consider one of the proposed designs for the LH2 tank of the Earth Departure Stage (EDS) for 
moon missions [23], in which the tank has the shape of a rounded cylinder with height Hq = 12 m 
and radius Rq — 2.5 m, containing 15 tons of LH2. The scale of the tank is similar to that used in 
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the S-IVB stage of the Saturn V rocket in the 1960’s and 70’s, and a brief comparison is, therefore, 
appropriate. 

Indeed, NASA’s most comprehensive experience with large-scale cryogenic storage tanks in 
orbital conditions goes back to the Apollo moon missions 1 [15, 16, 25] (see also reviews of other 
CFM experimental activities in [9, 14]). The third stage of the Saturn V rocket was propelled by 
LH2 and LOX, containing 19,800 kg and 88,800 kg, respectively (here and below the numbers are 
from the Apollo 17 mission [26,27]). About a quarter of the propellants (5,000 kg of LH2 and 
25,200 kg of LOX) was utilized for orbital insertion, and the remaining amount was used for the 
translunar injection burn. Once in orbit, the third stage spent about 3 hours in LEO. The tank 
insulation employed (polyurethane foam attached to the interior side of the tank wall) brought 
the LH2 boil-off amount down to about 1,000 kg, still an acceptable margin of under 10% of LH2 
available for the second burn. The flow of continuously vented hydrogen vapor (GH2) was used to 
provide enough thrust (on the order of 10 -5 g — 10 -4 g) to ensure that the propellants were settled 
at the bottoms of the tanks at all times during the orbital coast phase. The presence of small 
upward g-force ensured that the boil-off bubbles rose in the thin convective layer along the tank 
walls without entering the bulk liquid. Immediately prior to the second burn, the LH2 tank was 
rapidly pressurized by the stored heated helium gas (GHe), raising the tank pressure from about 
1.5 to 2 atm. This short-time pressurization, followed by firing of the ullage motors must have 
condensed the smaller vapor bubbles and made the larger bubbles move towards the ullage. The 
resulting LH2 liquid at the bottom of the tank should have, therefore, contained little or no bubbles, 
allowing to safely fire the engine. After that, the thrust of the engine would maintain the gas-free 
liquid at the LH2 intake, with screens adding an extra protection. 

While hugely successful in bringing man to the surface of the moon and back, this approach 
may not be applied to the missions currently under consideration. They key reason why the Apollo 
CFM approach worked for the lunar program was that the required in-orbit storage time for the 
cryogenic propellants was short enough, so it was possible to tolerate a large boil-off rate and, as a 
consequence, avoid microgravity conditions altogether during the time in orbit. Thus, the Apollo 
approach carefully avoided dealing with long-term CFM issues associated with microgravity. Any 
kind of large modern long-range mission, however, would require storing cryogens for extended 
periods of time. To achieve this, one would need to drastically reduce the amount of boil-off and 
work in micro-g or zero -g environments. New approaches are, therefore, needed to answer these 
emerging challenges. 

Modern multi-layer insulation (MLI) allows to dramatically reduce the boil-off rate compared 
to the Apollo missions. Let us assume that the tank is wrapped in an MLI blanket of 50 layers. 
Using the Lockheed correlation [28], we can estimate the heat flux through the MLI, given the outer 
environment temperature To ~ 240K, to be go — 7.4 x 10 -2 W/m 2 [29]. Taking for simplicity the 
tank area to be So = 2ttRoHo ~ 200m 2 , we obtain a lower bound of Qo = Qo&o — 15 W for the total 
heat flow into the tank, with the corresponding boil-off rate of at least 90 kg/month, or 0.6% of 
the propellant mass per month. Of course, these numbers must underestimate the actual heat flow, 
since they do not take into account heat leaks through various MLI penetrations by struts, feed 
lines, etc., as well as imperfections in the MLI itself. Let us also note that the heat flux depends 

x We note that liquid helium (LHe) has flown for extended periods of time on a number of scientific missions [13]. 
Let us point out, however, that in these missions LHe is chilled down to the superfluid state. This makes the case of 
LHe storage very different from all other cryogens, since, in contrast to all other cryogens, superfluid LHe has infinite 
heat conductivity [24], which prevents it from thermal stratification and nucleate boiling. 
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very sensitively on the outer environment temperature To. In the extreme case of To 400K, we 
find qo ~ 0.8 W/m 2 , an order of magnitude higher than the one computed previously, resulting in 
the heat flow of Qo — 160 W and an unacceptably high boil-off rate of 6% per month. 

Even with very efficient MLI insulation, the loss of propellant becomes prohibitive for extended 
missions. Therefore, active boil-off reduction techniques are necessary to improve retention of 
usable propellants. One idea that has been developed over recent years is to employ zero-boil-off 
technology (ZBOT) [3,30-32]. While ZBOT approach was demonstrated to be successful in the 
case of cryogens with higher boiling points, e.g. LOX, no cryocoolers enabling ZBOT yet exist that 
could operate at LH2 temperatures [33]. Another idea to further decrease the heat inflow into the 
LH2 tank is to use the concept of broad area cooling (BAC), whereby the tank is surrounded by a 
network of tubes carrying a coolant fluid [6, 12,31,33,34]. Circulating the fluid through the tubes 
with the subsequent heat removal by cryocoolers operating at higher temperature may then allow 
to significantly reduce heat penetration. Use of vapor-cooled shield (VCS) [6, 7, 13, 35] would be 
particularly efficient, since thermalizing GH2 with the outer layers of the MLI could increase the 
storage time up to a factor of 6 (see Secs. 3.1.4 and 3.1.5 for more details). We note, however, 
that strong localized heat leaks through MLI penetrations (the main subject of the present paper) 
provide one of the greatest CFM challenges for long-term cryogenic storage. 

3 Physics of long-term cryogenic storage 

Let us now perform some basic estimates for the sample 15-ton LH2 tank whose dimensions were 
introduced in Sec. 2. Since the precise parameters of the MLI performance admit significant 
variation, we will take an overall heat leak per unit area with an ample margin: qo = 0.4 W/m 2 , 
giving a total heat influx of 80W through the MLI (see also [6,29]). In addition, we will assume that 
another 40W of heat enters the tank via penetrations in the form of localized heat sources, giving 
the total incoming heat flow of Qo = 120W. We note that MLI penetrations, such as support struts, 
feed lines, etc., may provide the greatest challenge in the tank’s heat management. For example, 
taking the characteristic parameters of the orbiter support strut from the Space Shuttle external 
tank, which is a tubular structure of radius R ~ 20 cm, thickness d 5 mm and length L ~ 1 
m, with thermal conductivity k 7 W/(m-K) of Inconel 718 alloy at T ^ 100K [36], we find that 
the conductive heat leak into the tank can be estimated as Q stm t — ^m^RdT^/L — 11W. Note 
that this formula may significantly underestimate Q strut? since it does not take into account the 
additional heat entering the strut through its own thermal insulation. Similarly, for a titanium 
stmt with d ~ 1 cm and k 15 W/(m*K) [37] we find Qstrut — 45W, and for A1 2219 strut with 
d 1.5 cm and k ~ 70 W/(m-K) [38] we find a prohibitively high Qstrut — 320W. In view of 
the preceding considerations, however, the conductive heat leaks must not exceed several Watts 
per penetration in order for the tank to remain within the acceptable thermal budget. There is, 
therefore, a significant trade-off between the structural and thermal properties of the materials 
used, requiring strong materials with low thermal conductivity and a possible need for external 
penetration cooling [6]. 

To proceed, we need to specify the operating parameters for LH2 in the tank. We will assume 
that the tank is initially at pressure po = 1.6 atm, corresponding to the saturation temperature 
T s o = 22K of LH2 (parahydrogen), and that LH2 is subcooled to Tlq — 20. 3K, corresponding 
to saturation temperature at 1 atm. Here and everywhere below the definitions and values of 
the parameters of hydrogen used are are listed in Table 1 (for pressure po at saturation [39]). 
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Table 1. Physical parameters used in the estimates. 


Parameter 

Value 

Meaning 

D He 

5 x 10 -9 m 2 /s 

Diffusivity of helium in LH2 

Ho 

12 m 

Tank height 

P 

5 - 15 W 

Local heat leak power 

Q o 

120 W 

Total heat leak 

R o 

2.5 m 

Tank radius 

RgO 

1.8 m 

Ullage bubble radius 

Rh 2 

4,124 J/(kg-K) 

Gas constant of hydrogen 

RHe 

2,077 J/ (kg-K) 

Gas constant of helium 

So 

200 m 2 

Tank surface area 

T 0 

240K 

Exterior environment temperature 

Tl o 

20.3K 

Subcooling temperature 

T s o 

22K 

Saturation temperature of LH2 

Vo 

240 m 3 

Tank volume 

V 9 o 

24 m 3 

Ullage volume 

cl 

10,820 J/(kg-K) 

Specific heat of LH2 at constant pressure 

Cyj 

10 J/ (kg-K) 

Specific heat of aluminum 

9o 

9.81 m/s 2 

Earth’s acceleration of gravity 

9 

o - io-% 

Microgravity acceleration 

h 

1 cm 

Tank wall thickness 

Po 

1.6 atm 

Operating pressure 

9o 

0.4 W/m 2 

Heat flux through the MLI 

Ql 

4.35 x 10 5 J/(kg-K) 

Latent heat of LH2 vaporization 

Pl 

0.0192 K- 1 

Thermal expansion coefficient of LH2 

KL 

0.101 W/(m-K) 

Heat conductance of LH2 

K v 

0.019 W/(m-K) 

Heat conductance of GH2 

K"w 

20 - 200 W/(m-K) 

Heat conductance of aluminum 

PL 

1.16 x 10" 5 Pa-s 

Viscosity of LH2 

PL 

68.7 kg/m 3 

Density of LH2 

Pv 

2.07 kg/m 3 

Density of GH2 

Pw 

2,700 kg/m 3 

Density of aluminum 

CTL 

1.65 x 10" 3 N/m 

Surface tension of LH2 


Initially, a 10% by volume ullage space with volume V g o — 24 m 3 is present, pressurized by cold 
GHe. The required mass of GHe to produce the excess pressure of 0.6 atm is equal to M# e = 
(po — Patm)V g o/RHeTLo — 35 kg. We note that at large ullage volumes (as LH2 is lost due to 
boil-off or transfer from the tank), large amounts of cold helium gas are required to pressurize the 
tank. For example, when the ullage occupies 50% of the tank volume, one would need to supply 
175 kg of GHe, respectively. In practice, even greater amounts of helium may be required due to 
dissolution of GHe in LH2 on long storage timescales (see Sec. 3.1.6). The liquid is subcooled in 
order to avoid the presence of vapor bubbles in the bulk LH2, which are a potential hazard for 
engine restart, etc. 
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Let us briefly comment here on the nature of microgravity environment experienced by the 
cryogenic tank in LEO. For simplicity, we will consider the tank whose axis is oriented along its 
velocity vector in a circular orbit. The contributions to ^-forces can be separated into three parts. 
The first part has to do with the spatial non-uniformity of the earth’s gravitational field and the 
centrifugal force. The resulting effective g - force will be a linear function of the distance to earth, 
vanishing at the tank axial mid-plane. The apparent gravity g will point away from earth on the 
side of the tank farthermost from earth, and towards the earth at the side closest to earth. The 
maximum value of |g| can be shown to be given by |g| ~ 3go^o/^orbit — 10 _6 go> where go is the 
acceleration of gravity on earth. This gives the Bond number Bo = PlI^Rq/o'l — 2.5, indicating 
that capillary forces balance the apparent gravity forces in the case of the single ullage space. 
In this situation the equilibrium ullage shape should be approximately an oblate spheroid, whose 
center is located in the plane passing through the tank’s axis and is normal to the direction towards 
Earth. Note that more generally the apparent g- forces will be time-dependent on the time scale of 
the orbital period. From the balance of inertia and capillary forces, we find that the characteristic 
timescale for the ullage to settle is of the order of £ u n age \fpL^Jo~L - 15 mins. 

The second contribution to the ^-forces is the g-jitter, caused by the motion of the astronauts. 
These are impulses of acceleration with amplitude roughly of order 10 -4 go and duration 1 s in 
random directions (see e.g. the discussion in [40]). It provides the random vibrational background 
to the motion of fluid. Note that this contribution would be absent in an unmanned vehicle, or 
would be the only significant g - force on space missions beyond LEO. Also, both the g-jitter and the 
space-dependent apparent gravity may assist bubble coalescence. Finally, the third contribution to 
^-forces is due to the Coriolis force, which arises due to motion of the liquid inside the tank. 

3.1 Basic thermodynamics 

We begin by evaluating the heat budget of the tank and related issues. 


3.1.1 Time to saturation 


First, let us calculate the time needed for LH2 to come from the subcooled condition to the sat- 
uration temperature under an assumption of perfect mixing (e.g. by an active mixer inside the 
tank) and in the absence of any boiling and active cooling. This time is given by equating the total 
amount of heat that entered the tank to the increase in the LH2 heat content: 


= CLPL(F0 - y ^ ~ Tu>) * 26 days. 
Q o 


(i) 


On the other hand, in the absence of mixing, boiling, and any ^-forces the heat will only penetrate 
from the tank wall to the depth equal to the thermodiffusion length L = \J (clPl) of LH2 in 
time t. Then the same balance leads to 


t 


sat ,t hermodiffusive 


clPlkl(T s o - T LQ y 




16 days. 


(2) 


Note that, replacing T s o in Eq. (2) with 23K, one can see that in order to achieve a superheat of 
^1K at the tank wall (at which nucleate boiling normally occurs in LH2 under normal conditions on 
earth [41-43]), one would need to wait t ~ 40 days. Finally, for a tank in LEO, taking into account 



convective transport in microgravity with g = 10 -6 go, we obtain for the tank wall temperature T w 
an estimate T w —Tl q ~ 0.22K, which is based on the Nusselt number Nur 0 = q^Ro/ ( kl{T w — Tlq)) ~ 
46 and the Rayleigh number for these parameters Ra = g/3i(T w — Tlq)clp 2 l Rq/ — 2.8 x 10 7 , 
and we used the correlation Nur 0 = O.lSRa 1 / 3 [44, Eq. (9.31)]. Thus, the timescale on which the 
bulk of LH2 heats to the boiling point under the considered heat loads in the absence of any other 
sources and sinks of heat is about 1 month. Note that on such a long timescale heat conductance 
alone is sufficient to carry the heat into the tank interior. 


3.1.2 Time to complete evaporation 

On the other hand, assuming that the tank is maintained at constant pressure po> we can find the 
LH2 storage time by equation the total heat that entered the tank to the heat needed to vaporize 
all the LH2: 


w = = 22 months (3) 

Qo 

In particular, about 700 kg or ^5% by volume, of LH2 will be lost to boil-off in one month. This 
also means that cold GHe needs to be supplied at a rate of ^16 kg/month to enable maintaining 
LH2 at subcooled conditions. 


3.1.3 Self-pressurization and liquid heating due to mechanical work 

While it will take on the order of one month in the absence of any other heat sinks for the bulk of 
LH2 to be heated to the saturation temperature, local boiling may begin soon due to the localized 
heat sources through penetrations. With the power from those sources equal to Qo — qoSo = 40W 
and assuming that the boil-off bubbles remain attached to the hot spots on the tank walls, we can 
estimate the the rate of boil-off of GH2 as 


boil-off = — ° (T ° ° x - 8 kg/day. (4) 

Ql + c l (T s0 - T lo ) 

This results in an increase of the tank pressure at the rate of p = MGH2,boi\-offRH2T s o/VgO 0.3 
atm/day. 

Let us also point out that GH2 forming as a result of boil-off performs mechanical work against 
the liquid in order to expand. The mechanical work done on LH2 will be converted into heat 
through viscous dissipation in LH2, resulting in a further temperature increase in the liquid bulk. 
If the boil-off rate Mgh2, boil-off is given by Eq. (4), then the work done by the vapor on the liquid 
is 


w = PoMGg2 ’ boil - off ~ 7W. (5) 

Pv 

In other words, about 20% of the heat that goes into boil-off ends up heating the bulk liquid. To 
this, one should also add the amount of work supplied to the liquid through stirring. 
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3.1.4 Vapor-cooled shield and broad area cooling 

The use of VCS can significantly reduce the amount of boil-off by utilizing sensible heat of the 
boil-off vapor [6,7, 13,35]. If the GH2 is simply vented overboard, the amount of heat removed 
from the tank by the vapor will be Qboil-off — Ql^gH2^ where itigh2 is the mass flow of GH2 
through the vent. If, on the other hand, the vapor is allowed to thermalize with the outer surface 
of the MLI at temperature Tq 240K before being vented, the amount of heat removed equals 
Qvcs — (ql + cl(To — Tlo))i71gh2 (taking into account that the specific heat of GH2 does not vary 
significantly with temperature in the considered interval [39]). The ratio of the two is given by 


Qvcs cl (T 0 - T l o) 

— ~ 1 H ~ 6.5. 

Qboil-off QL 


(6) 


This means that the use of VCS may reduce boil-off rate for a fixed heat load by over a factor of 6 
by intercepting the heat entering into the tank. Further increase in the efficiency may be achieved 
by passing the warm GH2 through para-ortho converters [7] . 

Broad area cooling (BAC) is another promising concept that takes advantage of actively cooled 
GHe running through tubes within the MLI to capture the heat entering from the tank environment 
and allows to significantly reduce the heat flow per unit area through the MLI [33] (see also [34] 
for a related concept of active co-storage). One of the main difficulties in applying the BAC 
technology is efficient thermal bonding of the BAC tubing to the tank structures [6,31]. Recently, 
both the “tube-to-tank” and the “tube-to-shield” concepts, whereby the BAC tubing is bonded 
either directly to the tank wall or to an intermediate layer of the MLI, respectively, have been 
successfully demonstrated [12,31,33,45]. Let us point out that these concepts may also be used as 
part of a VCS in conjunction with TVS [35]. Note, however, that in the absence of boiling inside 
the tank the VCS/TVS system may not be able to intercept the incoming heat, since the heat 
entering the tank may raise the LH2 temperature locally near the tank walls without producing a 
pressure or bulk LH2 temperature rise. In this case, the TVS will not operate, and the propellant 
heating near the tank walls may lead to the potentially dangerous explosive nucleate boiling hazard 
(see Sec. 3.3.2 for more detail). We also note that in order for the tube-to-tank BAC design to 
work in an LH2 tank, the cold GHe must be circulated at the temperature of about 20K to avoid 
boiling of LH2 at the tank walls. Since no efficient cryocoolers currently exist operating at these 
temperatures, the tube-to-tank concept is not currently applicable to LH2 storage. 


3.1.5 Penetration cooling 

One of the most challenging problems in the tank heat management is to control the heat leaks 
from various MLI penetrations, such as the tank structural supports (struts) and propellant feed 
lines [6]. As was discussed in Sec. 3.1.4, using either VCS or BAC for passive or active cooling 
of the heat shield, respectively, one could all but eliminate the heat leaks through the MLI into 
the tank. In this case, the penetrations will provide most of the heat load to the propellant. Note 
that in contrast to the MLI, where radiative heat transfer dominates, it is not possible to adapt 
the active tube-to-shield BAC concept operating at intermediate temperatures (see Sec. 3.1.4) to 
intercept most of the heat coming through penetrations because of the dominant role of conductive 
heat transfer there. A natural idea (an extension of the VCS concept) is, therefore, to sacrifice 
some of the stored hydrogen to passively cool the penetrations at their points of contact with the 
tank walls and, in particular, to suppress possible boiling in those areas. The hydrogen is most 
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conveniently supplied by TVS, providing a regulatory feedback between the tank heat load and the 
amount of LH2 used for cooling. 

In view of the above discussion, to emphasize the physical phenomena under consideration, 
in this section we will assume that all heat entering into the tank comes from the penetrations. 
To make our considerations more quantitative, let us assume for simplicity that the tank contains 
Vgtrut = 8 identical penetration structures (struts). Consistently with our previous assumptions in 
Sec. 3, each strut is assumed to carry heat in the amount of Qstrut = Qo/^stmt = 15W conductively 
into the tank. Under our assumptions about the storage times for the tank, the maximum available 
LH 2 budget (not including the effect of VCS or BAC) that can be used to cool the penetrations is 

J stmt = = 3.3 x ltr 5 kg/s. (7) 

Ql 

This mass flow of cold GH2 and/or LH2 can be used, e.g., to locally cool the MLI penetrations in 
the shape of thick pipes. We consider a helical coil in perfect thermal contact with the strut surface, 
which is formed by winding a thin hydrogen-carrying tube around the strut (see Fig. 2). We now 
analyse efficiency of such a method of local cooling of the penetrations. In the following, we assume 
that the cooling tube has characteristic length L — 10 m. We also assume that TVS operates in 
the regime in which LH2 at po = 1.6 atm and Tl q = 20. 3K passes through a Joule-Thompson valve 
with pressure p\ — 0.2 atm at the other end. By considering an isenthalpic process, we then find 
that hydrogen will be exiting the valve in the form of a two-phase mixture at T\ 16K containing 
the mass fraction 77 0.9 of LH 2 [39]. We will consider the situation in which hydrogen from TVS 

is fully vaporized and is then thermalized with the LH2 propellant before being supplied into the 
cooling tube (see Fig. 2). In order to estimate the minimum tube radius i? t ube> consider that GH2 
is supplied in the form of gas at T — q. Using the Darcy- Weisbach equation [46,47] 


A _ fpu 2 L 
P 4 R tuhe 


(8) 


for the pressure drop A p across the tube, where / is the Darcy friction factor, p is the fluid density, 
u is the average fluid velocity, together with Haaland’s approximation to the solution of Colebrook’s 
equation for / [48]: 


/ = 


^1.8 logm 


/ trough V 11 ^9 

y 7.4i?tube / Rc 


-2 


(9) 


where d rcm gh is the tube surface roughness and Re is the flow Reynolds number, we can estimate 
the minimum inner tube radius i? t ube needed to supply the necessary LH2 flow J s trut given by Eq. 
(7). Taking d rcm gh — 1.5 pm [44] and A p — p\ — 0.2 atm, we find that the minimum tube radius 
is i?tube — 1-1 mm. Note that this estimate is quite insensitive to various assumptions on the 
parameters, such as the tube roughness or the pressure drop. 

Let us now estimate how effective these tubes will be in cooling the penetrations. Assume that 
the temperature of GH2 flowing inside the coil tube that winds around the hot strut changes from 
T v = Ti to T v — T 2 . By thermodynamic considerations, we must have 


Q strut — Cp(2~2 Ti)<Jstrut- 


( 10 ) 
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Figure 2. Cooling a penetration by GH2 from TVS. (a) The tank schematics and (b) the close-up 
of the cooled penetration. 


It then follows that the GH2 supplied at T\ = Tlq = 20. 3K will have the capacity to remove 
Qstrut = 15W from the strut, if it is heated to T 2 ~ 65K. Let us now estimate length of the tube 
which can remove this heat flow. For steady turbulent flow, the heat transfer coefficient /i tu be can 
be obtained from the Dittus-Boelter correlation [44] 

/Hube = ttA — N u tube , Nu tube ~ 0.023 x Re^ e Pr^ e , (11) 

^tube 

where Nu tll be> Retube and Pitube are the Nusselt, Reynolds and the Prandtl number associated 
with the gas flow in the tube, respectively. Note that the formula in Eq. (11) can be equivalently 
rewritten as 

(^)% r -V3, (12) 

— 9 /5 

i.e, at fixed J s trut the heat transfer coefficient is proportional to i? tub / e and, therefore, grows rapidly 
with decrease of Rtube- Considering a tube of radius i?tube = 2 mm, we find that for the mass flow 
given by Eq. (7) we have Re tu b e — 10 4 and Pr tu b e — 0.8, giving the Nusselt number Nu tu b e — 32 
and the heat transfer coefficient /i tu be — 130 W/(m 2 -K) (using the parameters for GH2 at T — Tlq 
and p = pi [49]). This corresponds to the power P loop ~ 47r 2 P tube P strut h tube (T s0 - T L0 ) = 1.8W 
taken away by the first loop of the tube, provided the strut at the base is at saturation temperature 
T s0 . Therefore, several loops are needed to take away the heat from the strut. Let us note that the 
pressure drop A p — pi ensures the flow J stru t i n a tube of length L ~ 200 m (see Eq. (8)). 

The analysis of heat exchange between the strut maintained at T = Tq at the warm end and in 
thermal contact with cold GH2-carrying tubes is rather involved and is presented in Appendix A. 
The results depend significantly on many factors, including the dimensions and the material of the 
penetration. Similarly, the ability of the tube to take away the heat from the penetration depends 
in a non-trivial way on these factors, as we demonstrate below. For concreteness, in the following 
we will assume (nominal regime) that a thermally insulated tubular strut of radius i? s trut = 10 cm, 
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Figure 3. Temperature distributions in the strut (blue line) and the tube (purple line) as a function 
of the distance along the strut to the tank wall, (a) Nominal regime; (b) The strut radius is 
increased to R s tru t = 15 cm; (c) The strut is shortened to T s trut = 30 cm; (d) The GH2 flow rate 
is decreased to J s trut = 2 x 10 -5 kg/s. 


thickness d s trut = 5 mm, length L stm t = 50 cm is made of pure titanium and take the value of 
^strut — 10 W/(m-K) for thermal conductivity, corresponding to the low temperature end of this 
material parameter [37]. Under these assumptions, the heat leak through the strut in the absence 
of vapor cooling is equal to Qstrut = 15W, consistent with the earlier discussion. We will also take 
the pitch of the helix dtube — 2 cm which is smaller than the strut radius and bigger than the tube 
diameter in order to ensure high surface area of thermal contact and absence of thermal shorts 
through the tubes. 

The solution of the governing equations for the temperature distribution in both the strut and 
the tube as a function of the distance from the tank along the strut under an assumption that the 
tube captures all the heat entering into the strut at the warm end is presented in Fig. 3(a). One can 
see that a temperature boundary layer develops near the tank wall, in which the strut temperature 
fails to follow cold GH2 temperature on the length scale /strut — 2.5 cm (see also the discussion in 
Appendix A). The maximum temperature at the strut base reaches 21. 6K, which is still below the 
saturation temperature T s o = 22K of LH2. Therefore, this rather crude analysis demonstrates that 
boiling at the attachment point of the strut may be suppressed by TVS-produced vapor cooling 
under suitable conditions. 

Let us now see how sensitive this conclusion is to variations of various parameters of the con- 
sidered cooling system. For example, consider the case in which the strut radius is increased by 
50%. The solution of the governing equations is then presented in Fig. 3(b). One can see that in 
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Figure 4. Cooling a penetration by LH2/GH2 mixture from TVS. (a) The tank schematics and (b) 
the details of the tube-to-tank scheme. 


this case the strut base temperature increases to T — 26. 5K, corresponding to a 4.5K superheat. 
Similarly, decreasing the strut length by 40% results in a strut base temperature T — 31K, a 9K 
superheat. Finally, a decrease of the GH2 supply rate by 40% results in the strut base temperature 
T = 29K, a 7K superheat. In all cases, the obtained superheat is big enough to initiate boiling. 
One should, therefore, carefully consider the heat conduction problem associated with penetrations 
in order to asses the feasibility of the proposed cooling strategy. One important factor to keep 
in mind is that, according to the model predictions, the efficiency of the proposed vapor cooling 
system decreases with an increase in the effective heat conductance coefficient of the strut. The 
latter must also include the effect of the additional heat conduction pathways introduced by the 
highly conducting thermal bonding material and the tubes themselves. It is not clear at this point 
whether this approach may always guarantee boiling suppression near the MLI penetration points. 
A more detailed and accurate analysis of the system is required. 

Before concluding this section, we briefly comment on the possibility to use LH2 from the Joule- 
Thompson valve directly rather than GH2 from the TVS heat exchanger to cool the penetrations 
(see Fig. 4). This can be accomplished, e.g., by the tube-to-tank configuration, carrying the 
two-phase LH2/GH2 mixture from the Joule-Thompson valve to the penetrations. Let us point 
out in the first place that the added cooling power from the latent heat of LH2 is in fact only 
a small fraction of the total cooling power of the GH2 stream. Indeed, for the same reason that 
the use of VCS may potentially increase the heat leak capacity of the LH2 tank by a factor of 6 
(see Sec. 3.1.4), the fraction of the latent heat of LH2 to the sensible heat of GH2, when heated 
from Tl o = 20. 3K to To = 240K is estimated to be qLJstrut/(cp(To — T^o)) — 0.2. Therefore, the 
advantage of using LH2 rather than GH2 to cool the penetrations is rather questionable. 

At the same time, in order to supply LH2 from the TVS output, the supply tubes must be 
sufficiently isolated to prevent boil-up and thermalization of the hydrogen flowing through the 
tubes with LH2 inside the tank. If not, LH2 will be transformed into GH2 due by the heat entering 
the tubes that carry the liquid from the Joule-Thompson valve towards the tubes wound around 
the strut (see Fig. 4). Indeed, the amount of heat entering a tube segment of length L may be 
estimated as K w Lh(TL o — Ti)/Z s , where k w is the tank wall thermal conductivity, h is the wall 
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thickness, and l s is the length of temperature equilibration due to convection. In the presence of 
natural convection alone in microgravity the value of l s may be estimated to be l s ~ 1 m (see 
Sec. 3.3.1 below for more detail), so for a highly conductive tank wall material with k w — 200 
W/(m-K) and h — 1 cm we get that LH2 in the supply tube will boil completely on the length 
L ~ rjqLJ stYll tl s / (rjK w h(TL o — Ti)) ~ 1.7 m. Therefore, bringing LH2 from TVS to a penetration 
without thermal insulation of the tubes on the way is problematic. On the other hand, suppose 
that by using suitable thermal insulation the heat entering the supply tube on the way to the 
penetration to be cooled is reduced sufficiently, so that LH2 at T\ — 16K is present at the first 
point of contact between the cooling tube and the penetration (see Fig. 2(b)). Then, since LH2 
in the tube is significantly cooler than the surrounding LH2 in the tank and the tank walls, the 
tube may draw a substantial amount of heat from the tank wall rather than the penetration. It 
is possible to adapt the analysis of Appendix C (leading to Eq. (18) below) to show that a 4K 
subcooling of a strut base of radius i? s trut = 10 cm would result in a power P ~ 17 W drawn from 
the wall, for the same tank wall parameters as above, in addition to P s trut = 15W coming from the 
stmt. This contradicts conservation of heat at the strut base. Therefore, the tube would not be 
able to cool the penetration, indicating that the tube temperature cannot be maintained at T = T±. 
This, in turn, implies the onset of film boiling in the cooling tube, which would greatly reduce its 
heat transfer coefficient to the values characteristic of the use of GH2 vapor. Thus, the advantage 
of using LH2 is once again greatly diminished. 


3.1.6 Helium dissolution hazard 

We point out that at Tlo = 20.3 K the solubility limit of GHe in LH2 is ^0.5% by weight [50]. 
Therefore, all LH2 in the tank is capable to absorb up to 75 kg of GHe, which is about 2 times 
more than the mass M# e ~ 35 kg of GHe in the ullage at the beginning. Let us estimate the 
time in which GHe may dissolve in LH2. Taking the diffusion coefficient of the dissolved helium 
DHe — 5 x 10 -9 m 2 /s (assumed to be of the same order as the available value for neon [51], see 
also [52]), assuming that the ullage has the shape of a spherical bubble of radius R g o = 1.8 m and 
estimating the diffusive flux at the ullage boundary to be P>HePHe,satl^He-, where pHe.sat — 0.005pL 
is the helium saturation density and Z# e = \jDn e t is the helium diffusion length, respectively, in 
LH2, we obtain that in the absence of active mixing the dissolved mass of GHe in time t is 

M-He, dissolved — / ^'^FlgO^HePHe,sat’) -^He, dissolved — 1-6 kg, t — 1 month. (13) 

Let us note that in the presence of active mixing the rate of GHe dissolution may be signifi- 
cantly higher. For example, if LH2 is circulated with an average axial velocity u^ik = 1 mm/s 
across the tank, then the Peclet number Pe = 2i? p o^&uZ/c / Dgje — 7 x 10 5 for the ullage bubble. 
Therefore, the dimensionless Sherwood number Sh = 0.65Pe 1//2 ~ 550 ( [53, Eq. (3.52)], assum- 
ing spherical ullage in a background flow with velocity u^ik), and the dissolution rate becomes 
MHe = ^RgoDH e pH e:Sa tSh 30 kg/month. Thus, the entire mass of GHe may dissolve in LH2 
in only one month, leading to the collapse of the ullage pressure. The dissolution rate may further 
increase due to ullage bubble motion. Therefore, on the long-term storage timescales one needs to 
evaluate the potential ullage collapse hazard due to GHe dissolution in LH2. One also needs to 
take into consideration the possible effect of dissolved non-condensible helium gas on the boiling 
characteristics [21,54,55]. 
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3.2 Pressure control 


When the sufficient level of superheat is reached at the tank walls and/or enough heat enters 
through penetrations, nucleate boiling will start. As vapor is created due to boil-off, the tank 
pressure rises, requiring venting in the absence of heat removal. As was shown in Sec. 3.1, with 
the considered heat leak the excess pressure po — p a tm will increase by a factor of two in just 2 days 
in the absence of venting. 

In microgravity, venting becomes and issue, since the location of the ullage space in the tank 
is not known and, therefore, it is not possible to guarantee that vapor, not liquid, is vented in 
the process. Note that venting LH2 caused the vehicle to tumble out of control during the AC- 
4 test flight in 1964 [8,9]. To circumvent this problem, the thermodynamic vent system (TVS) 
technology is proposed, in which the two-phase mixture of liquid and vapor may be expanded in 
a Joule-Thompson device and then passed through a heat exchanger to take away heat from the 
warmer fluid in the tank, see Fig. 5 for schematics [6, 10,31,56]. We note that ideally the work of 
TVS would result in the LH2 loss rate equal to that of a tank with the same heat budget venting the 
vapor directly overboard. In practice, increased losses are inevitable due to finite TVS efficiency. 

The basic physical principle of TVS operation is to control tank pressure by keeping the bulk 
LH2 at subcooled conditions. Therefore, successful operation of TVS relies crucially on its ability 
to efficiently move heat from the regions of nucleate boiling to the TVS heat exchanger inside the 
tank by utilizing fluid mixers, such as an axial jet or spray bars. We note that, on one hand, LH2 
circulation must be sufficiently strong in order for the heat through penetrations to be removed 
from vapor bubbles. On the other hand, if the circulation is too strong, then the ullage bubble may 
be distorted or fragmented, resulting in its possible capture by the TVS intake, see Fig. 6. 

3.2.1 TVS hazard induced by ullage motion 

In microgravity the ullage containing the saturated vapor and helium can drift to the TVS intake 
and become captured by it. This dangerous effect is unexplored. It may cause the TVS to cease 
normal operation. Bubbles containing helium may come out of the TVS, decreasing the helium 
mass in the main ullage. These effects may cause generation of many new bubbles with GHe and 
the formation of the multi-phase liquid-gas foam instead of the regular ullage. It maybe then be 
impossible to annihilate this foam by pre-start pressurization, since the collapsing bubbles will 
contain non-condensible helium gas. 

To estimate the critical LH2 circulation velocity Ub u ik, above which the ullage motion hazard 
becomes significant, we compute the flow Weber number We = 2pLul ulk R g o/aL. The flow will 
begin to affect the shape of the ullage bubble when the Weber number We > 1; for example, for 
U bulk — 0.5 cm/s we will have We ~ 4, indicating the onset of ullage motion [19]. 

3.2.2 Boiling near hot spots 

Let us now estimate the rate of heat removal from a bubble growing near a hot spot at the tank 
wall (Fig. 5(b) and Fig. 7). A detailed treatment of this question is difficult without specifics of 
the LH2 flow pattern, which depend on the tank design (see [32,57] for recent numerical studies). 
Therefore, for simplicity we will assume that LH2 flow is generated by the axial jet which creates 
a counterflow pattern (see Fig. 5), in which LH2 moves with average velocity Ub u ik in the region 
r < Rq — L, where r is the distance to the tank axis, and with average velocity u coun t er in the 
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|j\ Cold LH2 with temperature T L =20.3K 

' and pressure p=1.6 atm 

Interface temperature T S =22K 


GH2 
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Wall temperature T w > 22K 

Aluminum plate 



I I Hot spot with 

External heat flow q Q = 0.4 W/m 2 power P = 5W 


Figure 5. (a) A Tank/TVS configuration with an axial jet mixer, (b) Schematics of bubble growth 
on the tank wall. In (a), blue arrows indicate LH2 circulation; red arrows indicate heat influx 
through the MLI; on the left, a vapor bubble is attached to a hot spot (red) at the point of a 
localized heat leak through a penetration. 


opposite direction in a layer of thickness L next to the tank wall. Since at the wall the flow velocity 
is zero, in a laminar flow the magnitude of the LH2 velocity next to a bubble of radius R may 
be estimated as ui oc ~ 4 Ru coun t e r/ L, where we used linear interpolation from u — 0 at r — Rq to 
u = 2 u coun ter at r — Rq — L for the axial liquid velocity u. From mass conservation we find that 

4R(R 0 - L) 2 

Uloc ~ (2i? 0 - L)L 2 Ubulk 

Taking L~lm and the value of u^ik — 0.5 cm/s obtained in the preceding paragraph, one can 
see that a localized heat source with power P — 5W may produce a steady bubble of radius R~ 13 
cm. This follows from the balance of heat, taking into consideration that the Nusselt number 
Nu r = P/(2ttRkl(T s o — Tlq)) ~ 35 for the bubble can be correlated as Nur = 0.65Pe^ with the 
Peclet number Pe# = 2 ui oc RclPl/ kl — 3000 in this case (see [53, Eq. (3.52)]), assuming small 
contact angle and, hence, an almost spherical bubble, and ignoring the presence of the wall for the 
flow for simplicity). Similarly, a bulk velocity u^ik — 1 mm/s would result in a steady bubble of 
radius R = 20 cm. We note that even with this smaller value of the circulation speed u^ik one 
would need to supply LH2 through the axial jet at the rate Glh2 — k(Rq — L) 2 pLUi) U ik — 0.5 kg/s, 
which, for an axial jet with an orifice of 10 cm diameter would require the fluid velocity of about 
1 m/s. The corresponding circulation time for the whole tank would be 8 hours in this case. At 
the same time, for a given leak power the system will not be able to control the growth of bubbles 
attached to hot spots whose radii are smaller than the ones estimated above. 
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Figure 6 . An illustration of the TVS hazard induced by ullage motion (spray bar TVS schematics 
adapted from [56]). 


3.2.3 Formation of a multi-phase liquid-vapor foam 

The above estimates assumed that the bubble receiving power P from the heat leak remains at- 
tached to the hot spot at all times. In microgravity, bubbles will have a much reduced tendency to 
detach from the wall and rise towards the ullage [21,58-60]. Using the correlation of Fritz [17,61], 
we find that the bubble radius at departure in microgravity can be estimated at R d ~ 30 cm, 
assuming, e.g., a contact angle 7 ~ 20° (small contact angles are characteristic of LH2, a strongly 
wetting liquid [61,62]). Note that the Weber number for the bubble We# = 2pz J u‘f oc R/aL ~ 0.005 
is small in this case, and so the bubble may not be blown off easily from the hot spot by the flow. 
If the bubble departure radius Rd exceeds the steady bubble radius obtained above, the bubble will 
stop growing and assume its steady-state radius, and the heat will be transferred from the hot spot 
through the bubble to the liquid, as desired. 

The departure radius is proportional to the contact angle [17, 61] and becomes smaller than 
the steady-state bubble radius estimated in the preceding paragraph for contact angles 7 < 10 °. 
Thus, for these smaller contact angles bubbles will be leaving the area of the hot spot and entering 
into the bulk liquid. Other mechanisms of bubble departure, such as bubble coalescence, g-jitter, 
and the effect of the liquid flow [17,63,64], which are found to be important under reduced gravity 
conditions [58-60], may further reduce the bubble departure radius. As a result, bubbles may be 
injected into the liquid and start to move with the flow, reaching, in particular, the flow stagnation 
regions. Continuously arriving and collapsing bubbles will raise the temperature in the stagnation 
regions above the subcooling temperature T#o = 20 . 3K. As a result, the bubbles will stop collapsing 
and will accumulate in those regions. When bubbles arrive there, a complex multi-phase liquid- 
vapor foam-like mixture may form with temperature at saturation, T s 0 = 22K. The bubble colonies 
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Figure 7. The schematics of the tank wall in contact with a heat leak. A vapor bubble grows by 
partially absorbing heat from the hot wall. 


forming in this way may further thermally insulate the tank walls from the subcooled liquid and 
result in the formation of more bubbles through nucleate boiling, making TVS cooling ineffective. 

3.3 Dynamics of bubble growth and collapse 

We now discuss the dynamics of vapor bubbles inside LH2 in more detail. 

3.3.1 Onset of bubble nucleations at hot spots 

Let us begin by estimating the time needed for a single bubble to appear near a localized heat leak 
with power P = 5W. We assume that the heat enters the tank through a penetration connected 
to the exterior side of an aluminum LH2 tank. We note that in contrast to the Apollo design, in 
which thermal insulation was on the interior tank surface [16], the MLI has to be installed on the 
tank’s exterior surface, since it needs to operate in vacuum. Hence, in the absence of any special 
inner surface coating, the stored LH2 will be in contact with a highly conductive metal surface. 
At Tlq — 20. 3K, heat conductance k w of the tank wall lies in the range of 20 - 200 W/(m-K), 
depending on the composition of the aluminum alloy used [13,65]. Also note that bubble inception 
superheat will be higher than the one needed to maintain nucleate boiling and will vary depending 
on the tank wall finish [21,66-68]. 

Suppose first that LH2 is at the subcooled temperature Tlq when the heat leak is applied. 
Because of the much higher heat conductance of aluminum, this heat will first spread into the 
tank wall. If l w — -\Z^ w t/{c w Pw) is the thermodiffusion length of aluminum, then one can roughly 
estimate the temperature increase in the tank wall by equating the amount of heat Pt entering 
aluminum in time t to the heat content itl‘l 0 hc w p w {T — Tlq) of a cylindrical section of the tank wall 
with radius l w and thickness h. The timescale of temperature spreading may also be estimated by 
equating l w to tq\ 


t w = ~0.3s-3s, 


(15) 


where ro ~ 5 cm is the radius of the penetration. A more precise analysis (see Appendix B) gives 
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an extra logarithmic factor in the expression for the maximum temperature in the hot spot: 

T ~ T lo + — ^ — In (at/t w ), (16) 

4 nriK w 

where a ~ 6.1. Taking h — 1 cm, we find that for t — 3 t w the value of T varies in the range 21K 
- 26K, depending on the heat conductance k w . Therefore, three distinct scenarios are possible, 
depending on the values of k w and P in the absence of boiling. 

The first scenario is realized, if the heat leak power P is sufficiently high and the heat con- 
ductance k w is sufficiently low. For the considered value of P and k w — 20 W/(m-K), which is 
at the low end of the range of k w , by Eq. (16) we get T 24K already at t = t w = 3 s. This 
corresponds to superheat of 2K, already above the nominal IK superheat for the nucleate boiling 
onset in LH2 [41-43]. In this situation, nucleate boiling will start immediately, with all the heat 
going into a single vapor bubble. 

The second scenario is realized when the heat conductance tz w is sufficiently high. Then in 
the presence of losses through convection the heat from the leak spreads to distances up to (see 
Appendix C) 



where Nu# 0 is the Nusselt number associated with convection. The maximum temperature in the 
hot spot can be estimated as (see Appendix C) 

maxT ~ Tlo + 7 T~~r In (— ) , (18) 

27 TK w h V r 0 / 

where b ~ 1.85. For example, with k w — 200W/(m-K), corresponding to the high end of the range 
of k w , and Nur 0 ~ 46 estimated in Sec. 3.1 for free convection in microgravity, we obtain l s — 1 
m. About the same Nusselt number is also obtained for forced convection with average velocity 
U bulk — 1 mm/s considered in Sec. 3.1, using the correlation [44, Eq. (7.23)] 


Nu/? n = 0.332 x 


PLUbuikR ()\ 1 / 2 ( 


Then, according to Eq. (18), we get maxT ~ 21. 8K, so the temperature in a hot spot remains 
below the saturation temperature, and the heat is removed convectively from the leak, as desired. 
Note that in zero gravity and absence of any boiling or mixing the value of l s in Eq. (18) should 
be replaced with Rq, and the value of Tl q should be replaced with the spatially averaged tank 
temperature T w q. The latter will be increasing on the timescale t sa t 5 thermodiffusion (see Sec. 3.1), 
eventually leading to nucleate boiling. On the other hand, an addition of an active mixer next to 
a hot spot will further reduce the value of l s and, therefore, further suppress boiling. This may be 
a better strategy for mitigating the effect of heat leaks through MLI penetrations (compare with 
the strategies discussed in Sec. 3.1.5). 

The third scenario is realized when the wall heat conductance is low, and the heat leak power 
is also sufficiently low, provided that convective heat transfer is negligible. The latter takes place, 
e.g., in zero gravity and in the absence of active mixing. We consider this scenario in more detail 
in Sec. 3.3.2. 
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3.3.2 Explosive nucleate boiling hazard 

If the power P is not sufficient to initiate nucleate boiling quickly in time t w , then, according to 
Eq. (16) it may take an exponentially long time for the wall temperature to reach the required 
superheat, resulting in a long delay in the onset of nucleate boiling. Consider, for example, the 
case of Kj W = 20 W/(m*K) and P = 1W. Then, by Eq. (16) a temperature T ~ 24K corresponding 
to a 2K superheat will only be reached at t = 1.4 hours. In this time the heat will spread to 
l w ~ 2 m along the tank wall and II V ^Lt/ ( clPl) — 2.6 cm into LH2. When nucleation 
occurs, the heat stored in this superheated layer of LH2 will be used to convert a mass rriH2 into 
vapor, with rriH 2 — ^IwIlClPl{T — T 5 o)/^l = 1 kg for the considered parameters. The resulting 
vapor volume V — 0.5 m 3 will then be violently released into the tank in a short time. This is an 
example of the phenomenon of explosive nucleate boiling [67-69] that was ubiquitously observed in 
the microgravity boiling experiments on board the Space Shuttle Columbia [40]. We note that the 
accompanying pressure spike may present a potential hazard for the operation of the storage tank. 
Moreover, many vapor bubbles may be injected into LH2 as a consequence of explosive boiling, 
contributing to the formation of liquid- vapor foam. 

Let us note that in the absence of mixing the size of the hot spot is limited by the length scale 
(see Appendix D) 


L = ^h. (20) 

kl 

Assuming that k w 20 W/(m*K), as before, we find that L^2m. Correspondingly, the maximum 
temperature in the hot spot is limited by the expression given by Eq. (18) with l s replaced by L 
(see Appendix D). However, in this case the steady state superheat will extend into the liquid also 
to the length L, creating a much larger mass of superheated LH2, whose explosive boiling may lead 
to catastrophic consequences. 


3.3.3 Bubble growth over a hot spot 


Once a vapor bubble is nucleated, it will grow by drawing the heat from surrounding superheated 
liquid which, in turn, receives heat from hot spots on the tank wall. We note that for small heat 
fluxes considered here the dominant heat transfer mechanism will be transient heat conduction (for 
a recent discussion of different growth mechanisms, see [70]). 

Assuming that all the power from a localized heat leak is used to convert liquid into vapor in 
a single bubble, the bubble radius as a function of temperature, or, equivalently, the growth time 
for a bubble of a given radius, are given by 


( 3 Pt V /3 ^ _4irp v q L R 3 

\ Air f) v ([L J ’ growth 3 P 


( 21 ) 


With P — 5W, we then find 


R — 20 cm, £ growt h = 1.7 hour, R~ 2 cm, £ growt h = 6 seconds. (22) 

Note the difference in the dependence of the bubble radius Rout with the classical t 1 / 2 dependence 
obtained in [71,72] (see also [62,64]). This is due to the fact that the bubbles under consideration 
grow near a hot spot on a thin strongly conductive tank surface. Therefore, heat flux into the 
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bubble is mediated by the high heat conductance through the tank wall (for a recent discussion of 
the importance of heat conductance in the heater during nucleate boiling, see [73-75]; see also Sec. 
4.2). 

Note that far from the hot spots this formula remains valid, if one sets P — q^A, where A is 
the wall area per single bubble. The latter is valid in the case of high heat conductance of the tank 
wall, when on average each bubbles will be able to intercept all the heat entering through the area 
A of the tank wall. For example, for A = L 2 , with L — 10 cm, we find that the bubbles will reach 
the radius of R — 2 cm in t ~ 2 hours. When neighboring bubbles at the tank wall grow large 
enough to come in contact with each other, more complicated dynamics involving coalescence and 
detachment from the tank wall will occur. 


3.3.4 Bubble collapse and accumulation in the subcooled liquid 

As a result of several possible bubble departure mechanisms, vapor bubbles may detach from the 
tank wall and enter into the bulk liquid. Several scenarios are possible here, depending on the 
size of the departing bubbles, the level of microgravity, and the liquid flow rate. Larger bubbles 
may rise toward the area of zero gravity (the tank mid-plane for a tank in a circular LEO) under 
the action of buoyancy forces. The rise time may be estimated by balancing the buoyancy force 
PLg with viscous drag 4 tt^lRu, where u i?o Arise is the bubble velocity relative to the liquid 
(recall that Rq is the tank radius) [53], to obtain 


Aise — 


3/hlRq 

plqR 2 


~ 5 min, R — 2 cm, 


(23) 


which is decreasing with the increase of bubble radius R. These bubbles may then be swept towards 
the ullage bubble by the flow generated by the mixer and coalesce with it (see Fig. 5). The timescale 
of this process is given by 




ow 


H o 

Ubulk 


~ 3 hours, 


Ubulk = 1 mm/s. 


(24) 


Note that by Eq. (22) it takes much longer for a bubble to reach the ullage than to grow to the 
considered size. Hence bubbles may accumulate in the liquid and be carried along by the flow and 
deposited in its stagnation areas. 

Let us now estimate the bubble lifetime, assuming that upon departure it enters the region of 
the subcooled LH2. Assuming first that the heat escapes the bubble via steady conduction and 
equating the conductive heat flow 47tR 2 kl(T s o — Tlq)/R to the heat release rate 47 rR 2 (dR/dt)p v qL 
due to condensation, and then solving the obtained differential equation, we obtain the time for 
the bubble of radius R to collapse into the subcooled liquid: 

_ qlPvR 2 

icollapsel ~ 2kl{TsQ _ Tlq) 


(25) 


Note that this equation is asymptotically exact in the limit of vanishing subcooling, but underesti- 
mates t C oiiapsei f or larger subcoolings due to the condensation blocking effect [76]. Indeed, assuming 
that the condensation rate is dominated by the conduction through the thermal boundary layer of 
width l — yj kiAI{clPl) and, hence, one should replace the factor 1/R by l/l in the expression for 
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the heat flux (a more precise analysis, also giving the coefficient of proportionality is presented in 
Appendix E) 


, _ Kp 2 v q 2 L R 2 

tcollapse2 4 C L p L K L (T s0 - T L0 ) 2 ' 


(26) 


The blocking mechanism is not important when t co ii apse i < Aoiiapse2> which is the case in the 
examples that follow. 

From Eq. (25), one can see that a bubble of radius R — 20 cm will have a lifetime £ C oiiapsei 30 
hours (t growt h — 1.7 hour). For smaller bubbles with, say, R — 2 cm, we have Aoiiapsei — 20 mins 
(^growth — 6 sec). Note that these estimates agree very well with the results of direct numerical 
simulations of the full system of hydrodynamic equations describing the bubble collapse which are 
presented in Sec. 5. Once again, in the considered case the bubbles are produced faster than they 
are collapsing in the subcooled liquid, i.e. t gr0 wth Aoiiapse- Note, however, that as the bubbles 
collapse, they release heat into the liquid, so that the subsequent bubbles are exposed to a smaller 
degree of subcooling. Therefore, as more bubbles arrive, their lifetime may steadily increase and 
the liquid around them reach saturation, resulting in the formation of stable bubble colonies. These 
colonies then contribute to the formation of the multiphase liquid-vapor foam. 


3.3.5 Bubble collapse by pre-start pressurization 

The presence of vapor bubbles in LH2 is highly undesirable for the engine restart, since those 
bubbles may enter into engine feed lines and result in cavitation of the turbopumps [1]. A way to 
reduce or eliminate the vapor bubbles prior to engine restart is to pressurize the tank with GHe 
and, at the same time, apply a small settling thrust from ullage engines. Pressurization changes 
the LH2 saturation temperature relative to the bulk liquid temperature, making vapor bubbles to 
condense. 

Several issues arise in the course of pre-start pressurization which may result in an incomplete 
vapor bubble collapse, making the procedure inefficient. First, bubble condensation releases heat 
into the bulk liquid, increasing its temperature and potentially bringing it to the new saturation 
temperature and stopping further condensation. This is particularly relevant to the vapor bubble 
colonies. Consider, for example, the case in which the tank originally at pressure po = 1-6 atm 
is pressurized to a new pressure p — 2 atm. The corresponding new saturation temperature 
is T s — 23K. Now, suppose a bubble colony, which remains at the old saturation temperature 
T s o = 22K has vapor volume fraction /. Then the total amount of heat this foam-like multiphase 
fluid can absorb is 


Q = (1 - f)cLpLV(T s - T 5 o), 


(27) 


where V is the colony volume. This amount of heat, in turn, can condense only the volume 
V con d — Q/(qlPv) of the vapor. Comparing the value of V con d with the total vapor volume fV, we 
can see that by purely thermodynamic considerations all the bubbles will not be able to condense, 
if 


/> i + 


PvQl 


clPl{T s — T s0 ) 


-l 


0.45. 


(28) 


23 



In other words, it will not be thermodynamically possible to condense all the vapor, if the volume 
fraction of vapor exceeds a critical value given by Eq. (28). 

On the other hand, even if the volume fraction / of vapor is below the critical value, vapor 
condensation may not occur during the time interval of pressurization due to the general slowness 
of the condensation process and, in particular, due to the condensation blocking phenomenon for 
larger bubbles. One can once again use Eqs. (25) and (26) in these two regimes to estimate the 
collapse time for bubbles of different size, provided that T s o — Tlq is replaced with T s — T s q. Using 
these formulas, we now find that 

Wiapsei = 30 mins, R = 2 cm, Wlapsei = 50 hours, R = 20 cm. (29) 

It is clear that larger bubbles may not be eliminated by pressurization in a reasonable time. How- 
ever, by a settling acceleration g ~ 10 _4 go it is possible to move larger bubbles toward the ullage. 
Using Eq. (23) with Rq replaced by Ho, we find that bubbles of radii R > 5 mm will be able to 
move towards the ullage in time £ r i se = 5 mins. At the same time, for smaller bubbles the collapse 
time is bounded above by Aoiiapsei ~ 2 mins, assuming the best case scenario given by Eq. (25). 
This indicates that there are rather tight constraints for achieving the desired result from pre-start 
pressurization. Also note that during pre-start pressurization bubbles will be continuously gener- 
ated at the hot spots. Provided their departure radius is below 5 mm, these bubbles will not be 
eliminated at the moment of engine start. In addition, larger bubbles forming inside the engine’s 
start box, if any, will be trapped by the capillary screens of the liquid acquisition device (LAD) 
and will not be able to rise to the ullage. Finally, helium-filled bubbles forming as a result of a 
possible ullage capture by the TVS intake cannot be eliminated by pre-start pressurization. 


4 Mathematical models of bubble dynamics 

Here we present a discussion of more detailed mathematical models that can be used for the analysis 
of heat transfer mechanisms and bubble growth and interaction in microgravity under heat loads 
relevant for long-term cryogenic storage. 


4.1 Conductance-based models 


Consider a situation in which bubbles form due to boil-off at the tank wall in the vicinity of a 
heat leak in a quiescent saturated liquid (Fig 7). In the absence of vapor bubbles the heat from 
the localized heat leak will enter into the tank wall and rapidly spread along the wall. When the 
heat spreads sufficiently far from the source, it will then start entering into the liquid inside the 
tank. Mathematically, this can be described by the following equations (for simplicity, we neglect 
the temperature dependence of all the material parameters and assume that the liquid occupies a 
semi-infinite space next to the flat tank wall): 
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(30) 

z > 0, 

(31) 


where the first equation is the heat conductance equation in the wall, averaged over the wall 
thickness h (see [77] for a similar treatment), assuming the wall coincides with the plane z = 0, and 
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the second equation is the equation of heat conductance in the liquid. The two equations above 
are coupled via the boundary condition, namely, via continuity of temperature T at z = 0, and 
the derivative dT/dz is evaluated as z — ^ 0 + . The quantity q appearing in Eq. (31) represents the 
heat flux from a localized heat source. In writing the above equations we neglected the possible 
convective terms (including thermocapillary forces [54,55,58,59,78]), which are expected to be 
quite small under the considered heat loads. For example, in the case of a bubble of radius R = 2 
cm and heat penetration length l 2 i?, the associated Rayleigh number is 

R a = 0l(T ‘° - MSidll ~ 900, (32) 

f^L^L 

which is not sufficient for the onset of microconvection [79]. 


4.2 Bubble growth as function of bubble location 

Vapor bubbles play the role of heat sinks for the heat conductance problem in Eqs. (30) and (31), 
and will, therefore, screen the heat coming from the wall into the bulk liquid (for a recent review of 
single bubble heat transfer mechanisms, see [70]). To incorporate the bubbles, we need to introduce 
the boundary condition T = T 5 (p), where T s is the saturation temperature at the ambient pressure 
p, at the boundary of each bubble. As a consequence, heat will flow in or out of each bubble, making 
them grow or shrink, respectively. Under the considered conditions of very slow dynamics (on the 
time scale of hours or greater), the growth of a vapor bubble will be thermally-limited [71,72], 
with transient conduction dominating other possible heat transfer mechanisms [70]. For a spherical 
bubble (assuming zero contact angle for simplicity) this leads to 

dR f dT 

47 xp v q L R 2 — = Q, Q = K L — dA , (33) 

dt J dB du 


where dB denotes the boundary of the bubble, d/dv is the derivative in the direction of the outward 
normal to dB. Moreover, since the underlying conductance-limited bubble growth dynamics will 
be slow, it is reasonable to adopt a quasi steady-state approximation in solving Eqs. (30) and 
(31). This approximation will be valid when the time scale of thermal diffusion exceeds that of 
conductance-mediated growth. Estimating these timescales: 
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(34) 


one can see that the condition ^conductance ^ ^thermodiffusion is roughly satisfied for superheats under 
IK in LH2. 

Consider now an isolated bubble with the center ro which is at least distance d > R away from 
the wall. It is then possible to use the method of images [80] to express the solution of the obtained 
elliptic boundary value problem in the form of an infinite series: 


T(r) = T a (p) + 
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(38) 


where P is the mirror reflection with respect to the xy- plane. In view of the fact that q n +2/qn < 
R/d < 1, the series in Eq. (35) converges exponentially fast. Furthermore, by Gauss theorem, the 
total heat flow Q into the bubble is given explicitly by 


Q = - 2 



<72n(r) d 2 r. 


(39) 


This expression may then be substituted into Eq. (33), which gives an equation for the bubble 
radius in terms of the flux q from the wall to the liquid. This flux, in turn, can be found by solving 
the integral equation obtained by plugging the wall temperature, obtained from Eq. (35), into Eq. 
(30) with the time derivative set to zero. 

To get a sense of the phenomena described by the equations obtained above, let us consider a 
situation, in which a large bubble is suspended in the liquid at a distance d from the wall, such 
that R <C d « L, where L is defined in (20), directly over a point heat source q(x,y) = P5(x)5(y). 
Then it is possible to retain only the first term in the series in Eq. (35) . This allows one to estimate 
Q as follows 


Q ~ 



PRd 

"L2"’ 


(40) 


where we took into account that q ~ P/L 2 and that only the region of size of order d will give a 
significant contribution to the integral in (39). According to Eqs. (33) and (20), this implies that 
the bubble will grow on the time scale t ~ p v qLR 2 h^h 2 / (PdK 2 L ). Taking R d ~ 1 cm and P 
W, we can estimate t ~ 10 hours for k w = 20 W/(m-K). On the other hand, if d R, the entire flux 
P will be absorbed by the bubble, giving a much faster time scale of growth t = p v qLR? / (3P) 
(i.e. Eq. (21)), and for R = 1 cm we now have t ~ 4 s. Thus, the timescales associated with the 
bubble dynamics depend sensitively on the bubble location. 


4.3 Interaction of bubbles 

From the arguments leading to Eq. (40), it is clear that many bubbles may need to be created to 
take away the heat entering the tank from a localized heat leak. In principle, it is possible to use the 
approach of Sec. 4.2 to incorporate the effect of multiple bubbles on heat transfer. However, with 
the number of bubbles rapidly increasing, the problem of calculating the series representation of the 
temperature distribution quickly becomes intractable. Instead, a homogenization approach which 
reduces the bubble-filled liquid to an effective medium is advantageous (for a general reference, see 
e.g. [81]). Note that a related class of problems was recently treated by homogenization techniques 
in [82]. 
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Figure 8. Multiple bubbles attached to the wall (a) or inside the liquid (b). 


There are two main situations that need to be considered (see Fig. 8). In the first, it is assumed 
that only bubbles attached to the wall are present (Fig. 8(a)). In this case each bubble will work 
as a sink absorbing a heat flow Qr — Crk w H(T — T s )/\n(l/R), where l R is the characteristic 
distance between bubbles and Cr is a constant that needs to be obtained from the solution of 
the homogenization cell problem. Introducing the distribution function / 2 (i?, r,£) which gives the 
number of bubbles of radius between R and R + dR in the wall area element d 2 r around r at time 
t, we then find that the averaged heat flow into the bubble per unit wall area is 

poo /nr 

Q(r,t) = K 2 (r,t)(T(r)-T s ), K 2 (r,t) = K w h — f- f 2 (R,r,t)dR , (41) 

J o i n (V R) 

where K 2 is the effective (homogenized) heat transfer coefficient in two dimensions. Substitut- 
ing this formula into Eq. (30), we now obtain the respective homogenized equation for the wall 
temperature 


dT _ 
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If all the bubbles have the same radius i?, and the spatial density of bubbles is constant equal to 
C 2 = 1/Z 2 , then the coefficient K 2 becomes 


_ 2c2K w hCR 
2 |ln(c 2 i? 2 )|' 

From this equation, one can obtain the screening length A 2 y/hK w /K 2 l ln 1 / 2 (//i?). This will 
be the length scale of spreading of heat from a heat leak in the presence of bubbles attached to the 
wall. 

In the second scenario one needs to consider bubbles suspended in the liquid. The homog- 
enization problem in this situation (with fixed bubbles) was solved in Ref. [83]. The resulting 
homogenized version of Eq. (31) becomes 

dT (d 2 T d 2 T d 2 T\ 

™ St =KL {w + W + w)- h ' 3{t ■ t)(T - T - ] ' (44) 
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where, introducing the distribution function fs(R, r, t) of bubbles of radius between R and R + dR 
in the volume element d 3 r, we get 

poo 

K 3 (r,t) = 4 tck l Rf 3 (R,r,t)dR. (45) 

Jo 

Once again, if c 3 = l/l 3 is the constant density of bubbles inside the liquid, we find 


K 3 = 47 tk l Rc 3 , 


(46) 


and the respective screening length is A3 V ^jKl-l^l/R. 

One can further use the expressions obtained above to find the effective heat transfer coefficient 
K 2 appearing in Eq. (42). Assuming that T varies on a much longer spatial scale than A3, we can 
reduce Eq. (44) to a one-dimensional boundary value problem (also taking advantage of the quasi 
steady-state approximation). As a result, the solution for T (homogenized) in the liquid may be 
written as 


/ poo \ 

T(x, y, z)~T s + (T(x, y, 0) - T s )e~ z / X 3 , A 3 = 4vr / Rf 3 (R)dR 
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(47) 


where we assumed for simplicity that the bubble distribution is constant in space and time. Plugging 
this expression into Eq. (30), we obtain Eq. (42) with 


K 2 = 


KL 

As’ 


(48) 


giving the screening length in the plane A2 ~ yj hX^t^w / — VTA3. The problem needs to be 
studied further when both mechanisms of heat transfer act simultaneously. 


4.4 Ostwald ripening 

Let us now discuss a possibility of Ostwald ripening in the suspension of bubbles inside the liquid or 
a collection of bubbles attached to the wall [66] . In microgravity, an intriguing new mechanism may 
result in Lifshitz-Slyozov-type dynamics of vapor bubbles, even in the absence of mass transport 
between bubbles. Instead, the mechanism can be mediated via temperature, as we demonstrate 
below. Due to surface tension, the vapor pressure inside a bubble of radius R is equal to p — 
Po + 2 gl/R [84]. At the same time, in order to achieve evaporation-condensation equilibrium, the 
vapor pressure must be equal to the saturation pressure at temperature on the bubble wall. This 
leads to a peculiar feedback mechanism: The temperature on the surface of a smaller bubble with 
radius R\ is higher than that on the surfaces of a larger bubble with radius R 2 : 

Ti > T 2 if R, < R 2 , Ti , 2 = T s + 2 Ts Z L , (49) 

qLp v R\,2 

where we used the Clapeyron-Clausius relation [84]. As a consequence, heat will flow from the 
smaller bubble to the larger one, resulting in condensation in the smaller bubble and evaporation 
in the larger bubble under suitable conditions. It is possible to show that the resulting problem 
is mathematically equivalent to the one describing Ostwald ripening during solute precipitation in 
a supersaturated solution and can, therefore, be described within the Lifshitz-Slyozov theory [85]. 
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Estimating the time scale of ripening by balancing the amount of heat An R 2 p v qndR needed to 
increase the bubble radius by dR with the inflow of heat 47 tR 2 klT s gl/ { qLp v R 2 )dt in time dt, we 
obtain 


^ripeningl — 


klolTs ' 


(50) 


and the usual Lifshitz-Slyozov R ~ t l ' ? ' growth law. For LH2, we find that this process is quite 
slow in the case of bubbles suspended in the liquid, with t r ipeningi ^2.5 days for bubbles of average 
radius R ~ 1 mm. At the same time, if the bubbles are attached to the highly conducting tank wall, 
then the heat conductance between the bubbles is enhanced, with the heat flow into the bubble of 
order 27 rRhK w T s crL/(qLp v R 2 )dt, resulting in an estimate 




ripening 2 hn w O L T s ' 


(51) 


and an unusual dependence R ~ t 1 / 4 . Once again, using h — 1 cm and k w = 200 W/(m*K) for 
highly conductive tank wall material, we obtain t ~ 20 days for R — 2 cm. Therefore, the latter 
ripening mechanism may play an important role in selecting the bubble sizes. 


4.5 Motion of bubbles suspended in the liquid 

Let us now discuss the motion of bubbles under the action of microgravity. Let g(r,t) denote the 
apparent acceleration of gravity in the tank (which includes forces of inertia). Note that both the 
direction and the magnitude of the vector g may vary in space and time. Since the magnitude of 
g is very small, the movement of bubbles will be dominated by viscous drag. An isolated bubble 
away from the walls will, therefore, move with velocity [53, Eq. (3.15)] 


v = u — 


PlR 2 
3 Ml 


g, 


(52) 


where u is the background liquid velocity, and we neglected the vapor density compared to pl- 
Note that when |g| ~ 10 -5 m/s 2 and R ~ 1 cm, the bubble would move with speed |v| ^ 2 mm/s 
in LH2, comparable to the characteristic fluid circulation velocities considered in Sec. 3.2.2. In the 
case of many hydrodynamically non-interacting bubbles, one can use this equation for v, together 
with Eq. (33), to write the continuity equation for the bubble density fs(R,r,t): 
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(53) 


The heat flow Q , in turn, depends on R and the temperature distribution in the liquid, obtained, 
e.g. by solving the homogenized Eq. (44). In writing Eq. (53) we also introduced an effective 
diffusion constant D depending on R which possibly arises due to turbulence or vibrations (the 
g-jitter is not included in g). Finally, note that the bubbles, in turn, will exert a body force onto 
the liquid and may, therefore, affect the flow. 
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4.6 Bubble creation and departure from the wall 

Equation (53) needs to be supplemented with boundary conditions at the tank wall. Let us note 
that because of microgravity, the bubbles that form at the surface by heterogeneous nucleations 
will not easily detach from the wall. One source of such detachments are vibrations. At the same 
time, another important mechanism was recently identified in the experiments conducted in the 
low gravity environment of the NASA’s KC-135 aircraft [58-60]. In this mechanism, two bubbles 
have to be sufficiently close to each other, so that at a certain moment their coalescence occurs. 
The energy released in this process may be sufficient to break the resulting large bubble free from 
the wall [60]. Taking these two effects into account, we can write the following boundary condition 
for fs in Eq. (53): 
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(54) 


where k\ and k 2 are the corresponding single and double bubble detachment rates, and f 2 is the 
distribution function of the bubbles attached to the tank wall. The latter, in turn, is governed by 
its own balance equation which contains growth terms and sink terms due to bubble departure: 


t + M (rap h)=M(R-Ro)(l-A„f MR, r, t )dR ) 

r oo /♦oo 

-fci(i?)/ 2 (iW)- / / k 2 (R, i? 7 , i? // )/ 2 (i? / , r,t)f 2 (R'\ r, t)dR r dR ,r , (55) 

J 0 J 0 

where J is the nucleation rate per unit area, which depends on the local temperature T(r,£), Aq 
is the area per nucleation site, and Rq is the radius of newly nucleated bubbles. The first term in 
the right-hand side of Eq. (55) states that all empty nucleation sites will be activated with rate J. 


5 Computational analysis 


We now illustrate how the basic physics analysis discussed in the preceding sections can be supple- 
mented by high-fidelity computational studies of the underlying fluid dynamics equations in order 
to obtain accurate information about the bubble dynamics. For definiteness, we consider the prob- 
lem of vapor bubble collapse in subcooled LH2 upon bubble departure from a localized hot spot. 
To simplify matters, we assume that upon departure a spherical bubble is carried along with the 
flow of the subcooled liquid, where it shrinks uniformly by vapor condensation. 

Let us begin with the governing equations for the liquid phase in the vicinity of the bubble. 
The motion of the liquid phase can be described by incompressible Navier-Stokes equation, written 
under the assumption of spherical symmetry [79]: 


du 2 u f du du\ dp f d 2 u 2 du 2 u\ 

dr r ’ \dt dr ) dr \ dr 2 r dr r 2 J 5 


r > R, 


(56) 


where r is the radial coordinate, R = R(t) is the bubble radius, u = u(r,t) is the radial liquid 
velocity for r > R, and p — p(r, t) is the pressure inside the liquid. The liquid velocity u vanishes 
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at infinity. Similarly, the thermal diffusion equation in the liquid has the form [79] 


'dT dT\ 

CLPL l — + U— ) = K L 


dt 


dr J 


d 2 T 2 dT 
dr 2 r dr 


r > R, 


( 57 ) 


where T — T(r, t) is the liquid temperature for r > R, and we neglected the heat generation due 
to viscous dissipation. From the conservation of mass, the corresponding boundary condition for 
the velocity u at r — R are 


j(t) = p L (u(R + 0,t) - R(t)), 


(58) 


where j is the condensation mass flux at the bubble interface, and here and below the dot denotes 
time derivative. 

Since the considered bubble dynamics is slow compared to the timescale of sound propagation 
in the vapor phase, with very good accuracy the pressure p v inside the bubble is independent of 
space (a similar situation is realized in models of gaseous combustion [86]) and can be obtained 
from the condition of mechanical equilibrium at the interface: 


Pv(t) =p(R,t ) - ~ j 2 (t)(p v l (R ,t) -Pl 1 ), (59) 

where the second term accounts for viscous forces [53], the third term is due to a pressure jump 
due to surface tension [53], and the fourth term is the recoil force due to evaporation/condensation 
flux [75,87]. With this assumption in mind, let us now turn to the governing equations for the gas 
phase. The conservation of vapor mass is described by [79] 


dp v d(p v u ) 2 p v u 

dt dr r 


r < R, 


(60) 


where p v — p v {r,t) is the vapor mass density and u = u(r,t) is the radial vapor velocity for r < R. 
Furthermore, neglecting the contributions of the gas velocity and treating the vapor as an ideal 
gas, we can write the energy conservation in the form [79] 


d(^vPyT) djppPj/up) 2cpPyuT_ _ f <9 2 r 2 dT\ 

dt dr r Kv \ dr 2 r dr J ’ 


r < R, 


(61) 


where T = T(r, t ) is the vapor temperature for r < R and c v is the specific heat of vapor at constant 
volume. Note that the ideal gas equation of state used in Eq. (61): 


Pv(t) = p v (r , t)R v T{r , t), r < R, 


(62) 


relates the vapor density p v and temperature T at every point inside the bubble. The continuity 
of mass and energy at the interface read 

j it) = p v (R,t)(u(R — 0,t) — R), (63) 


and 


dT(R — 0, t) dT(R + 0,t) . 

QLj(t) = Kv ~ kl a . T(R — 0, t) = T(R + 0, t), 


dr 


dr 


(64) 
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respectively. Finally, to close the obtained system of equations, we relate the condensation flux 
j with the conditions at the interface via the Hertz-Knudsen relation [88,89], using a simplified 
equation for saturated vapor pressure p s as a function of the temperature at the liquid-vapor 
interface [90]: 


= Pv{t) -p s (T(R,t)) 
v /27ri2 w T(i2, t) ’ 


Pa(T) =p c 



(65) 


where R v — Rh 2 , Pc = 1-315 x 10 6 Pa and T c = 33. 2K are the critical pressure and temperature, 
respectively, and A = 5 for hydrogen. 

The continuity equation in both the gas and the liquid phase allows one to solve explicitly for 
the fluid velocity u at every point in space. For the liquid phase, which is incompressible, the 
solution is well-known: 


/ x a(£) 

u{r,t) = 


r > R, 


(66) 


where a(t) is some function of time only. Then, integrating Eq. (56) over r, with the help of Eq. 
(66) we obtain 


, a a 

PyR) t) PO PL ( ^ — 2^4 


Also, substituting Eq. (66) into Eq. (57), we have 

( dT a dT 


CLPL (-m + ^fr' =KL 


d 2 T 2 dT 
dr 2 ^ r dr 


r > R. 


(67) 


(68) 


In the gas phase, we can use the assumption that the vapor pressure p v is space-independent 
to obtain a similar equation for u. Indeed, from Eqs. (60) and (62), we find that 


du 1 dT u dT p v 2 u 

dr = T~dt + T~dr ~p^~~' 


r < R. 


Substituting this expression into Eq. (61) and using Eq. (62) yields 


(dT dT\ 

C ^{m +U ^) =K 


d 2 T 2 dT\ . 

v l Tpy H -X - + P v ■ 

dr z r dr J 


(69) 


(70) 


Now, let us use Eqs. (69) and (62) in Eq. (60), multiply the obtained equation by r 2 and integrate 
over r. After some algebra, we obtain a formula for u : 


" (r, ‘ )= s4®( 3<7 “ iK ^y r<R • 

where 7 = c p /c v . Combining Eq. (71) with (70) then yields 

c pPv dT k v fdT\ 2 rp v dT_ f d 2 T 2dT\ . 
R^T~dt + Y\dr) 3(7 - \)T~dr ^ ^ ydr* + r~dr) +Pv 


(71) 


(72) 
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Figure 9. The dynamics of bubble collapse in the subcooled liquid obtained by solving the governing 
fluid dynamics equations: (a) the bubble radius R and (b) the condensation flux j, as functions of 
time. 


We now substitute the obtained solutions for u into the boundary conditions at r — R. First of 
all, from Eqs. (58) and (66), we find that 


£ = 4 -— • 

R PL 

On the other hand, from Eqs. (62), (63) and (71) we obtain 


Pv = 


R 


7 Pv 

PL 


(7 - 1 )cpT(R,t) 


7 Pv a , ( dT(R,t ) 

+ (7 - l)«i 


i? 2 


dr 


(73) 


(74) 


where j is related to other variables via Eq. (65). Now, using Eq. (59) together with Eqs. (66) 
and (67), we also have 


R , 

H \Pv~Po~ 


2 R 3 PL 


4 Plcl 2u l 


R 3 
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+ f 


(7 - 1 )cpT(R,t) 


IPv 


Pl 


(75) 


This equation closes the system of equations describing T{r,t), R(t), a(t ) and p v {t ). Thus, in order 
to obtain the full solution to the bubble dynamics, one needs to solve Eqs. (64), (65), (68), (72) - 
(75). For the problem of bubble collapse, the following initial conditions may be used: 


-R(O) = Rq, 

Pv(0)=po, a( 0 ) = 0 , 

( 76 ) 

T(r, 0 ) = T s q, r < Rq, 

rp( rp | ( 2 "s 0 77 . 0 ) 7^0 ^ p, 

2 {r, 0 ) = l L o + , r > R 0 . 

nr* 

( 77 ) 


The results of the numerical solution of the above equations under a simplifying assumption 
that p v (t) = p s (T(R(t),t)), justified in the considered situation [76], and with Rq = 2 cm are 
presented in Fig. 9. As expected, the bubble shrinks to zero radius in finite time due to vapor 
condensation on the liquid-vapor interface. The heat released from the condensation is carried 
away into the subcooled liquid, and the temperature inside the bubble stays essentially constant, 
see Fig. 10. One can also see that the bubble collapse time is in very good agreement with the 
estimate obtained in Sec. 3.3.4. 
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Figure 10. Temperature (a) and velocity (b) distributions in and around the collapsing bubble at 
times t — 0, 300, 600, 900 s. 

6 Summary and conclusions 

In conclusion, we have presented an overview of the physics-related issues associated with large- 
scale, long-term storage of LH2 in microgravity. A special difficulty associated with LH2 storage, as 
opposed to other cryogenic liquids, such as liquid oxygen or liquid methane, is its very low boiling 
point (^20K). Currently, no efficient cryocoolers exist that can operate at LH2 temperatures. This 
presently puts the promising ideas of ZBOT technologies out of reach in the case of LH2. For the 
same reason, other active cooling technologies, such as BAC, are also limited to the “tube-to-shield” 
concept, which can only partially address the problem of LH2 heat load reduction. Importantly, in 
the case of LH2 storage heat leaks through MLI penetration may become dominant challenges for 
enabling long-term storage in microgravity. 

Below we first list the main physics-related issues of LH2 storage and then discuss potential 
challenges identified by our analysis. 

In microgravity, the location, structure, and shape of the ullage space is generally unknown. 
The dominant forces determining the ullage bubble shape are capillary forces. The location of the 
ullage bubble is determined by a combination of the ^-forces and the forces exerted on the ullage 
bubble by the liquid stirring flow. The ullage bubble may drift inside the tank due to the generally 
time- varying character of the ^-forces. 

Since the ullage bubble position is unknown and may be variable, boil-off vapor cannot be vented 
directly overboard. A combination of TVS and a fluid mixer may be used to control the boil-off 
pressure by maintaining the liquid in a subcooled state. This requires an addition of non-condensible 
cold pressurant gas to avoid ullage bubble collapse. Further heat load reduction strategies, both 
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active (BAC applied to an intermediate layer of the MLI) and passive (VCS utilizing the boil-off 
vapor from TVS to cool the MLI and penetrations) may also be necessary. 

The stirring flow velocity of the liquid inside the tank is subject to competing constraints. 
Stirring has to be sufficiently strong to ensure efficient heat transfer from the tank walls through 
the liquid into the TVS heat exchanger. Stirring has to be sufficiently weak to avoid deforming or 
breaking the ullage bubble. Excessive stirring may result in the injection of unwanted energy into 
the liquid. 

Ullage motion driven by a strong stirring flow may result in a TVS hazard, whereby the ullage 
bubble may be captured by the TVS intake and block normal operation of TVS. As a consequence, 
TVS may fail to control the tank pressure, or bubbles filled with non-condensible pressurant gas 
may come out of the TVS recirculation outlet. 

This implementation faces a number of challenges related to fundamental physics. 

1. On the basis of correlations derived under earth gravity, it is conceivable that both natu- 
ral convection in microgravity and stirring may keep the tank wall temperature below the 
saturation point away from the possible temperature hot-spots due to localized heat leaks 
through MLI penetrations. Efficient heat removal from the penetrations may require highly 
conductive tank wall material of sufficient thickness and vapor cooling. On-orbit experiments 
and further physics-based modeling is necessary to establish new reliable correlations which 
are needed to design a proper strategy for tank temperature control. 

2. Unpredictable changes of the flow conditions may result in overheating of the liquid near the 
tank wall. For high wall conductivities this may lead to the explosive boiling hazard, whereby 
a large mass of liquid becomes superheated in a rather large area around a temperature hot- 
spot on the tank wall or near flow stagnation points after a long delay time. The ensuing 
fast nucleate boiling may create a pressure jump and a strong liquid flow affecting the ullage 
bubble. 

3. Single large vapor bubbles attached to the tank walls or many small bubbles injected into 
the liquid may form at the points of strong localized heat leaks. The choice of the scenario 
depends sensitively on the mechanism of bubble departure from the tank wall during nucleate 
boiling near hot-spots. Key factors determining this mechanism are the liquid contact angle, 
the tank surface properties, microgravity magnitude and direction, g-jitter, Ostwald ripening 
and bubble merging. 

4. The bubbles inside the liquid may drift with the stirring flow and accumulate in its stagnation 
areas, resulting in the formation of complex multi-phase liquid-vapor foam-like structures. 
The foam may further thermally insulate parts of the tank wall from the subcooled liquid, 
increasing the boiling rate in those regions. The associated pressure increase cannot be well 
controlled by TVS. 

5. Boil-off bubbles may not be easy to remove by an overpressurization, since the collapsing 
bubbles release heat into the surrounding liquid and raise its temperature to saturation. In 
the presence of TVS hazard by the ullage motion, overpressurization may be ineffective in 
the presence of non-condensible pressurant gas inside the vapor bubbles. 

6. The multi-phase liquid-vapor foam may grow at the expense of shrinking of the ullage bubble, 
when the vapor mass is transferred from the ullage into the foam. This phenomenon is 
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complicated by the pressurant gas dissolution hazard, which is further amplified by strong 
stirring. The dissolved non-condensible pressurant gas may enter into the vapor bubbles of 
the foam. As a result, these bubbles may be difficult or impossible to eliminate by applying 
overpressurization. 
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Appendix A 


Analysis of heat exchange between the cooling tube and the strut 


Here we derive a set of approximate differential equations that describe the distribution of 
temperature in a tubular strut of material with thermal conductivity /%trut 5 radius R s t ru t, thickness 
^stmt 5 and length L s tru t, with a helical coil made of a GH2-carrying tube wound around (see Fig. 2). 
The pitch of the tube helix is denoted by d t ube and is assumed to be much smaller than i? stm t- Let 
x > 0 be the coordinate along the strut and let x — 0 correspond to the point of the strut attachment 
to the tank wall. Consider the heat balance in the stmt on the interval [xq — ^dtube^o + l^tube]- 
Assuming that the stmt temperature T(x ) varies appreciably only in the x-direction, to a first 
approximation the steady heat balance reads 

— ^ 0 _ l _ 2 ^t;ube 

27ri?strut^strut^strut = 47T ^Rstrut^tube^tube (^(^o) ^u(^o))? (A78) 

X x=x Q -\d tVLhe 

where T v (x ) is the vapor temperature at position x along the stmt in the tube, /i tu be is the heat 
transfer coefficient obtained in Eq. (12), and we assumed that the tube is in perfect thermal 
contact with the stmt. Using Eq. (11) to express /i tu be i n terms of the tube Nusselt number and 
approximating the jump in the derivative of the temperature above by dt^cPT / dx 2 evaluated at 
x = xq, after some algebra we obtain the following differential equation for T(x): 


‘'strut 


d 2 T 

dx 2 


= T-T„ 


/strut — 


' ^strut ^strut <Aube 
7T /V^Nutube 


with boundary conditions 


dT 

dx 


= 0 , 

x =0 



x — -^strut 


(A79) 


(A80) 


where we took into account that all the heat from the stmt is captured by the tube and thus 
does not enter into the tank at x = 0 , and that the stmt is thermalized with the environment at 
x = L stru t* Note that the parameter /strut has the dimension of length and characterizes the length 
scale at which the stmt temperature would equilibrate to the tube temperature, if the latter were 
kept constant. For the parameters of Sec. 3.1.5, we find Z s trut — 2.5 cm. 

On the other hand, the heat balance for GH2 in the interval [xq — ^dtube^o + ^dtube] reads 


c pJs trutTy (x) 


X — Xq-\- 2 ^tube 
X=Xq 2 ^tube 


4tT -R s t ru t-Rtube^ , tube(^~'(^o) ^u(^o))- 


(A81) 


Once again, approximating the jump in T v (x) by d tu b e dT v /dx at x = xo, we can write this as a 
differential equation 


t 7 dT v 
CpJ strut ^tube 


4tT -R s t ru t-f^tube^ , tube(2~'(^o) ^u(xq)). 


(A82) 
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Then using Eq. (A78) and integrating the obtained equation, we get 

rri rri l7 _ 27ri? s t ru t d s trut ^strut / a QQ\ 

— -L LO \ G , 5 t'v — j > (^AooJ 

CpJgtrut 

where we took into account that GH2 enters with T v = Tl q at rr = 0. Here the parameter l v 
also has a dimension of length and characterizes the length scale on which the vapor temperature 
approaches that of the strut, if the latter follows that of the vapor (when cooling is efficient). For 
the parameters of Sec. 3.1.5, we find l v 10 cm. 

Finally, substituting Eq. (A83) into (A79), we obtain a single linear equation 

d 2 T dT 

^ tTut ~dx 2 + ^ v ~dx ~ T ^ L0 = (A84) 

Together with the boundary conditions in Eq. (A80) this simple boundary value problem admits 
an exact solution, which is too cumbersome to be written down as an explicit formula. We omit 
the final answer, which is plotted in Fig. 3 for various parameters. 


Appendix B 

Analysis of heat spreading in the tank wall 


Consider a flat infinite layer of aluminum of thickness /i, to which heat power P is injected 
uniformly into a disk of radius ro > h at one of the surfaces. Averaging the heat conductance 
equation written in cylindrical coordinates over the film thickness, we obtain an equation 


d ( dT\ 


dT ~ , 

c wp W gt ~ r Qr p Qr J + 


7 ITq/i 


ff(ro-r), 


(B85) 


where H(x) is the Heaviside function. The solution T(r, t) of this equation with initial condition 
T(r, 0) = Tl q can be written in terms of the two-dimensional heat kernel (see e.g. [91]). From that 
solution, we find that the maximum temperature, which is equal to the value of T(0,£) is 


Tmax(t) 


Tlo + 
Tlo + 
Tlo + 




pt dl r fro 

/ 7 — 77/ rdr exp(-c w p w r 2 /(4K w (t -t'))) 

Jo t ~ t Jo 

CwPwTq \ 1 


47 ic w p w K w hrl 

P 


AttIiKo, 


In 


4 




c w pwr Q 

1 - e ^ wt ) + c w p w r 0 T 0, 


4 K w t J 


e 7 l c w p w r£ 


+ 0{t 




(B86) 


where T(a,z) is the incomplete Gamma-function and 7 0.5772 is the Euler constant [92]. Ap- 

proximating the solution in the tank wall by Eq. (B86) and evaluating the coefficient inside the 
logarithm numerically, we obtain the estimate in Eq. (16). 


Appendix C 
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Effect of heat loss into the liquid on the heat conduction in the 
tank wall 


Consider an infinite layer of aluminum of thickness h injected with a heat leak P uniformly 
through a disk of radius ro > /i, as in Appendix B, in contact with a layer of liquid of thickness 
Rq undergoing free convection. If Nur 0 is the Nusselt number associated with convection and 
klNur 0 /Rq is the corresponding heat transfer coefficient, then the equation of heat conductance 
in the layer, averaged over the layer thickness (see e.g. [77, Appendix] for a similar treatment) and 
written in polar coordinates is 


dT _ 
CwPw ~di ~ 


d 


r dr 


dT\ 

r lfr ) 


(T — Tip) 
hR o 


+ 


irr^h 


H(r 0 - r). 


(C87) 


In contrast to Eq. (B85), this equation admits steady state solutions as t —$■ oo, due to the screening 
effect of the heat loss, characterized by the screening length given by Eq. (17) (similar equations 
arise in the studies of Debye- Hiickel theory of electrolytes, in which l s is the Debye radius [84]). 
Indeed, the stationary solution of Eq. (C87) is [93, Eq. (3.39)]: 


T(r) = T lo + 


P 

irr^Kwh, 


Is - rol s Ki(r 0 /l s )Io(r/l s ), r < r 0 , 
rol s Ii(ro/l s )Ko(r/l s ), r>r 0 , 


(C88) 


where I n (x) and K n (x) is the modified Bessel functions of order n of the first and second kind, 
respectively [92], and so the solution decays exponentially at lengths of order l s . The maximum of 
T(r) is achieved at r = 0, and we have 


maxT(r) = T L0 + 


P 

2l TK w h 


In 


2 exp (| - 7 ) l s 


r 0 


r 2 1 

+ 0 (^ln^ 

Ij r 0 


(C89) 


where 7 ~ 0.5772 is the Euler constant [92]. Evaluating the numerical factor inside the logarithm 
and dropping higher order terms, we obtain Eq. (18). 


Appendix D 


Spreading of heat from the wall into the liquid 


Consider a flat infinite layer of aluminum at z — 0 covered by a semi-infinite layer of LH2 at 
z > 0. The steady heat equation for aluminum averaged over the layer thickness reads 


'q2 T q2 T \ Kl g T 

x '" ’ dx 2 dy 2 / h dz 


+ ^k H ( r °-^ T ^) =0 ’ 


z =+ 0 ^ r l h 


(D90) 


where the last term is as before and the boundary term comes from the continuity of heat flux 
from LH2 into aluminum. The corresponding heat conductance problem for LH2 is given by the 
following boundary value problem: 

d 2 T d 2 T d 2 T 

+ TTT z > 0 

ox z oy z oz z 
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and continuity of T is assumed down to z = 0. Applying Fourier Transform in the x and y variables 
to Eq. (D91): 


/ +oo /»+ oo 

/ e iqix+iq2V (r(a , 5 z ) _ Ti0 ) dx dy , 

-oo J — oo 

where q = (gi, ^ 2 ), and solving the obtained equation, we get 

f q (z) = T q (0)e-^, 


q = Jqi+qI 


(D92) 


(D93) 


In particular, the Fourier Transform of the Dirichlet-to-Neumann map associated with (D91) is 


dT n 


dz 


= -?r q ( 0 ). 


(D94) 


z =+0 


Plugging this expression into the Fourier-transformed Eq. (D90) and using the expression of [93, Eq. 
(B6)] for the transform of the characteristic function of a disk, after some algebra we obtain 


r q ( 0 ) = 


2PJi(qr 0 ) 
r 0 K L q 2 (Lq + 1)’ 


(D95) 


where J n (x) is the Bessel function of the first kind of order n, and we used Eq. (20). Inverting 
the Fourier transform in Eq. (D95) and performing the integrations, after some algebra we find 
ma xT(i, y, z ) = T(0, 0, 0), with 


maxT = Tjj } + 


2 n L r 0 


*(?)-*( 


2 L 
7rr 0 


(D96) 


where H\{z) and Yi(z) are the Struve function and the modified Bessel function of the second kind, 
respectively [92]. Finally, expanding the expression in Eq. (D96) in a series in ro, we find that 


maxT = T lo + 


27 TKonh 


In 


a 

^0 


2exp( 5 - 7 )L'| +0 /r| ln L 
} \L 2 r 0 


(D97) 


which coincides with Eq. (C89) with l s replaced with L, to the leading order. 

Appendix E 


Bubble collapse in the subcooled liquid via conduction through 
thermal boundary layer 

Consider a bubble of radius i?(£), whose interface is at saturation temperature T s o suspended in 
a subcooled liquid with temperature Tl o- In the thermal boundary layer approximation the heat 
conduction near the interface may be approximated by a one-dimensional initial boundary value 
problem for the liquid temperature T(r, t) 

dT d 2 T 

c lPl~q ~ = K L~r^ 2 1 T(R,t) = T s0 , T(oo, t) = T lo , T(r, 0) = T L0 . (E98) 
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The solution of this problem is 


T(r, t ) = T lo + ( T s0 — T lo ) erf 


( [cLpiJr - -R ) 2 
4 K L t 


(E99) 


where erf (z) is the error function [92]. From this formula, we find that the heat flux through the 
interface is 


dT(R,t ) / clplkl(T s o Tlo ) 2 

9(f) = = V 


(E100) 


The substituting this expression to the heat balance equation 4irR 2 dRp v qL = — 47r R 2 q(t)dt, whose 
solution is given by Eq. (26). 
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